Approximating Steady States in Equilibrium and Nonequilibrium Condensates 
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We obtain approximations for the time-independent Gross-Pitaevskii (GP) and complex GP equation in two 
and three spatial dimensions by generalizing the divergence-free WKB method. The results include an explicit 
expression of a uniformly valid approximation for the condensate density of an ultracold Bose gas confined in a 
harmonic trap that extends into the classically forbidden region. This provides an accurate approximation of the 
condensate density that includes healing effects at leading order that are missing in the widely adopted Thomas- 
Fermi approximation. The results presented herein allow us to formulate useful approximations to a range of 
experimental systems including the equilibrium properties of a finite temperature Bose gas and the steady-state 
properties of a 2D nonequilibrium condensate. Comparisons between our asymptotic and numerical results for 
the conservative and forced-dissipative forms of the GP equations as applied to these systems show excellent 
agreement between the two sets of solutions thereby illustrating the accuracy of these approximations. 
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PACS numbers: 31.15.xg, 03.75.Hh, 47.37.+q 

The complex Gross-Pitaevskii (cGPE) equation, also 
known as the cubic Ginzburg-Landau equation, or Nonlin- 
ear Schrodinger equation (NLS) arises in many branches of 
physics. It has successfully been used to model phenomena 
such as nonlinear (optical) waves, second order phase transi- 
tions, superconductivity, superfluidity, and Bose-Einstein con- 
densation of atomic gases as well as quasiparticle excitations. 
In the context of superfluidity, the GP HI equation has served 
as an excellent model for atomic gases while the cGPE has 
faithfully reproduced a number of key phenomena observed 
in experiments on nonequilibrium condensates of quasiparti- 
cle excitations. A key feature of this equation is that it de- 
scribes phenomena dominated by different physical processes 
that lie on either side of a nonlinear turning point. The non- 
linear turning point is governed by a second Painleve tran- 
scendent, which is a canonical equation arising in all of the 
contexts mentioned above, and more generally in systems that 
are nonlinear generalizations of an underlying linear prob- 
lem governed by a second order differential equation. It is, 
therefore, no surprise to see that it also arises in the nonlinear 
Landau-Zener problem fl. Yet, unlike Airy's equation, which 
governs the classical turning points of the linear Schrodinger 
equation, a uniformly valid solution or approximation for this 
equation has remained a formidable challenge. 

Focusing on the problem of Bose-Einstein condensates for 
ultracold gases, we note that the most commonly used method 
for determining the steady state solutions of the (GP) equa- 
tion is based on the Thomas-Fermi approximation. However, 
in general, it is not possible to extend the Thomas-Fermi (TF) 
profile uniformly into the classically forbidden region be- 
cause of the need to solve the second Painleve transcendent 
J51 01 ■ The difficulty in obtaining uniformly valid approxi- 
mations has meant that the TF approximation has been used 
in many circumstances where it is clearly invalid. A particu- 
larly important example arises in determining the equilibrium 
properties of a Bose gas at finite temperature. To determine 
the equilibrium properties of a weakly interacting Bose gas 
within a confining potential, Nikuni and Griffin [5] invoked 
a TF profile for the condensate density, and a WKB (equiv- 



alently a local density approximation (LDA)) for the thermal 
excitations. This in turn produced an unrealistic cusp shaped 
distribution of the thermal cloud density at the edge of the con- 
densate. Since the thermal cloud attains its maximum value at 
the edge of the condensate, the error introduced by the TF ap- 
proximation introduces significant errors in the computation 
of the thermal cloud density. 

Aside from its relevance to experiments, knowledge of the 
equilibrium properties of the system is also especially impor- 
tant in finite temperature models of Bose gases. A common 
feature of many of these works (e.g. c-field methods J^j]) is to 
model the macroscopically occupied coherent part of the sys- 
tem using a classical field that is coupled to an incoherent part 
of the system that is made up of higher energy scarcely oc- 
cupied modes. An energy cut-off must then be specified that 
determines which subset of the system is modeled as a clas- 
sical field. This cut-off is typically determined from the equi- 
librium properties of the system y7j] requiring a solution of the 
Bogoliubov-de-Gennes (BdG) equations in the Bogoliubov or 
Popov approximation. Being able to specify the cut-off in a 
simpler way through useful analytical approximations would 
be particularly useful for such finite temperature models. 

In addition to the above problems, there has been a surge 
of activity in recent years in the properties of nonequilib- 
rium condensates. These can include exciton-polariton Ji| 
or magnon ]9[] condensates where a condensate is created 
through coherent pumping to balance the dissipative processes 
that exist in such systems. The action of these nonconserva- 
tive effects can significantly alter the form of the condensate 
from the TF profile as was illustrated by Keeling and Berloff 
ifloll . Given these recent developments, useful approximations 
that go beyond the TF approximation are clearly needed. 

In ll 111 , a method was proposed that resolves the diver- 



gences that arise around the turning points in the classical 
WKB methods. In contrast to other approaches 012I1 . the di- 
vergence free WKB method also provides a uniformly valid 
solution for the ground state wavefunction of the NLS equa- 
tion. However, the divergence free WKB method has found 
limited applications partly because it is restricted to ID sys- 
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tems. Since Bose-Einstein condensation is also studied in sys- 
tems of higher dimensions, we will begin by extending the re- 
sults of the divergence free WKB to higher dimensions. For 
generality, we consider the cGPE given by 

itldrtfr = -L-Vfy+Vaaffr + gMfy + HS-DW, (1) 
2m 

where S and T> denote pumping and dissipation respectively. 

To begin, we will focus on the T — Bose gas where 
S = T> — 0. In this case, the equation has two constants of mo- 
tion corresponding to a fixed total energy which in 3D is given 
by the Hamiltonian H = f (Ti 2 /2m\V^\ 2 + V ext \if/\ 2 + § W 4 )d 3 x 
and the total number of particles N = J\if/\ 2 d 3 x. The param- 
eter g = Anash 2 jm, where a s is the s-wave scattering length 
for the interatomic interaction potential. In this work, we will 
be concerned with a harmonic trapping potential of the form 
Vext = j(<^ 2 x 2 + a>yy 2 + a> 2 z 2 )- Following 11311 . we will non- 
dimensionalise Eq. (Q~|) using the energy scale fuoho and the 
length scale R = ah (15Noa s /aho)^ 5 , where No are the num- 
ber of particles in the condensate (note No+ N for finite tem- 
perature models based on the c-field approximation), a/ 10 = 
yl%l(moJho), and the average oscillator frequency is defined in 
terms of the three oscillator frequencies as w/,„ = (oijWjW;) 1 ' 3 . 
Letting i/r -> VM/re"'' 1 ' 2 leads to 

id t if/ = -—V 2 i//+V en i// + yW 2 ^-fi^. 



(2) 



where V ext (x) = \{AlJ+%y 2 +Ay),y = 4nN a s a A JR 2 , e 2 = 
(aho/R) i s a small coefficient, is the chemical potential, and 
the trap anisotropy is given by A x = cj x / ajho, A y = iD y lu>ho, A z = 
a> z /u>i 10 . Despite the nonlinearity in our equation, we proceed 
by seeking variable separable steady-state solutions. Moti- 
vated by the divergence free WKB for ID systems, we express 
the wavefunction as t//(x,y,z) = exp[(ip(x) + §(y) + (f>(z))/e]. 
Substituting into Eq. (0, we obtain 



o = T W 



-&" + 4>")-Uf' 2 + ^' 2 + <P' 2 ) + V, 



(3) 



Differentiating the above equation with respect to x, we obtain 



0: 



f ff 

■<P<P 



+ ^YBL +2y '£ e l {ip+ ^)le_ (4) 
dx £ 



After eliminating <p" with Eq. ((3), we have 



e-ip 



■ = if' [e(§" + <[>") + <p' 2 + §' 2 + <p' 2 + 2(^1- Vext)] + e 



3VW 
dx 



If the trapping potential is of the form V ext (x,y,z) = Vi(x) + 
Vi(y) + V?,(z) (as for a general harmonic trap), the above equa- 
tion is variable separable. It then follows that we can write 



eV" 



<P 



(5) 



where n c = ji + C. From the boundary conditions d x tf/ = d y tf/ = 
d-ifr = at x = (0,0,0) we have (p'(0) = #'(0) = f(0) = and 



so the constant C = e(§"(0) + (p"(0))/2. By similar reasoning, 
one can arrive at analogous equations for &'(y) and <p'(z) re- 
ducing our 3D problem to a solution of three uncoupled equa- 
tions. We can now proceed as in the divergence-free WKB 
to obtain solutions for <p' by neglecting the terms of 0(e 2 ) on 
the left-hand side of Eq. © whilst retaining terms of 0(1) and 
0(e) on the right-hand side. As discussed in 111 111 , the resulting 



cubic algebraic equation of the form, <p + aip + btp' + c — 
and similarly for #',0', has three roots one of which corre- 
sponds to a uniformly valid solution of the cubic NLS equa- 
tion. If a,b, and c are real, this root is given by 



</t=(A + B)--, A = -sgnOR) 



Q/A, (A*0) 
(A = 0) 



Q 



\R\ + 



-3b 



yjR 2 -Q 3 



1/3 



2a i -9ab + 27c 
54 ' 



The solution we have obtained for \j/(x,y,z) = 
exp [f(p'dx + J§'dy + ffidz], can therefore be evalu- 
ated in terms of a quadrature of the integrals in the exponent. 
This leading order solution includes quantum mechanical 
effects (e.g. healing layers) arising from the kinetic energy 
terms. Higher order corrections can now be obtained to this 
zeroth order solution by expanding <p' in powers of e such 
that <p' = zZT = o e <Pn- At the next order, we obtain 



(6) 



Now for a harmonic trap, it turns out that <p' Q can be integrated 
explicitly. Using integration by parts together with Eq. Q to 
express x in terms of <p' , we obtain 



X 



r ip' Q dx = xip' Q - J x?d<p' =x(p' -^ In \<p' Q \ + i [e 2 + fy c tp Q 2 

+*p'o] 1/2 + y { ln + fo) + (e 2 + ^' 2 + 4v' 4 ) l/2 }} 

-^arctanh{(6 2 +4 f i cl p' 2 )/[e(e 2 + 8^ 2 +4 V ' 4 ) l/2 ]\ - w (0). 



As far as we are aware, this explicit expression for <pt pro- 
vides the first uniformly valid explicit approximation for the 
condensate density that includes the effects of the healing lay- 
ers near the edges of the condensate. Moreover, this solu- 
tion applies to 2D or 3D harmonic traps with different os- 
cillator frequencies in each spatial direction allowing cigar 
and pancake shaped condensates to be obtained. At next or- 
der, the analytical solution requires the evaluation of a simple 
quadrature for the integral J Q ip'^dx' . The solution, presented 
above, determines (p upto some constant. To specify the con- 
stant, we evaluate the value of ^(0) from Eq. (f3]i It follows 
that <p(0) + 0(0) + m - f In { ^'(O^WCO)]**' } . We note 
that explicit dependence on the parameter y appears only in 
the value of ^(0). However, its effect is also contained in 
the chemical potential fi through the normalization condition 
J \i//\ 2 d 3 x=l. 

Figure QJ presents a comparison of the condensate density 
for a spherically symmetric harmonic trap against numeri- 
cal simulations for the ground state of the GP equation for 



3 



H = 23.05. By specifying // to correspond to the value in our 
numerical simulations, our asymptotic solutions will approx- 
imately satisfy the normalization condition on the wavefunc- 
tion. The numerical simulations were performed using La- 
guerre polynomials as the basis functions. We see that even 
at leading order, the results are in excellent agreement with 
a fidelity of 99.957%. The small discrepancies that arise can 
be reduced if we include the contribution ip' x arising from the 
next order in our expansion as seen in the inset. The numeri- 
cal and analytical results almost fully coincide at this order of 
approximation with a fidelity of 99.996%. The second inset 
shows the functions tp' , and tp' x illustrating that higher order 
corrections are localized around the nonlinear turning point. 

We have also extended the solution to a rotating conden- 
sate with a vortex located at the center of the trap by seek- 
ing radially symmetric solution. Now following 11411 . we 
first transform the radial NLS equation into a form suitable 
for a WKB approximation by setting tfr(r,z) = U(x,z) where 
r = e x . We then seek a variable separable solution of the form 
ifj(r,z) = U(x,z) = exp([ip(x) + <f)(z)] + isd) where 9 is the polar 
angle. Following a similar procedure as discussed above, we 
can reduce the problem to two equations for <p' and <p'. The 
equation governing the radial wavefunction for the case s — + 1 
is now given by 



+ e V ' 2 + (2fie 2x - e 4v - e V + ee 



4.v 



— e = e 



2 if 



-<p 



In order to truncate the above equation, we perform a local 
analysis around the origin in which the terms containing the 
exponentials become exponentially small. Proceeding with 
this simplification and neglecting the terms on the right hand 
side, we obtain <ffi + eip' 2 - e 2 ip' - e 3 =0. This equation has 
two solutions given by w' = ±e which can also be expressed in 
terms of <p'(r) = <p'(x)/r = ±e/r. The solution near the origin is 
therefore given by U(x) = U(r) = e f ±e/r ' dr ' = e ±elnr = r ±€ . On 
physical grounds, we obtain U (r) = r e as expected. Moreover, 
if we now substitute this solution back into the right hand- 
side of the full equation given above, we find that these terms 
vanish near the origin. This justifies why we can neglect these 
0(e 2 ) terms whilst retaining the 0(e 2 ) and <9(e 3 ) on the left 
hand side which allows us to satisfy the boundary conditions 
for the vortex at the origin. We, therefore, compute tp' by 
neglecting only the right-hand side of the above equation. 

Proceeding as before, we can then obtain an explicit expres- 
sion for <p and (p. Whilst <p has the same form as the expression 
given above, <p(r) is now given by 

rw' 1 

$>(r) = r\n(r)ip' + — - ~ eMW + el) 



-I in 

2 



firijp' - V(^') 4 + (M 2 ~ 2e 2 )(rip') 2 + e 4 



arctan 
2 



rip' - e 

fi[(rip') 2 + e 2 ] 



rip 



2ey/(rip') 4 + (n 2 -2£ 2 )(rip') 2 
e l + (rip') 2 + 



e 4 ) 



2e 2 ){rip') 2 + e 4 



where ??(•) denotes the real part. The above expression for 
ip(r) is known upto some arbitrary constant that sets the overall 
normalization condition for the wavefunction. In this exam- 
ple, we have computed this constant such that the maximum 
value of \if/\ corresponds to the maximum value obtained nu- 
merically. We have included a comparison of the computed 
and analytical profiles in Fig.QJ. As before, we see excellent 
agreement even in this case of a rotating condensate with a 
quantized vortex located at the center of the trap. 

Having illustrated our approximation of the condensate at 
T — 0, we now extend the results to a Bose gas at finite tem- 
perature. We recall that the BdG system in the Bogoliubov 



approximation can be written as 115 



16] 



— V z (D(x) + (Vext + y«o)<i>(x) 



-8 l 2 

—V ui(x) + (y ext + 2yno-p)ui(x) - ynov/(x) = 
-5 2 , 

—V v ; (x) + (Vext + 2y«o-/')v;(x)-ynoM;(x) = -Em(x), 

where fi = jjticJi w I 'kgT ', «o(X) = M)|0(x)| 2 is the density of 
the condensate, and iit(x) = Ei?t0^i(l M i( x )l 2 + l l, !'( x )l 2 ) is me 
density of non-condensed particles. These equations are 
solved subject to the conditions N - Nq + 2iV0^i where Af; = 
{exp [(£, -/*)] - l}" 1 for i + 0. We have non-dimensionalised 
our equation for the condensate <J> as before for the T — case. 
The equations for the Bogoliubov modes, ui and v, have been 
non-dimensionalised using the length scale R and the energy 
scale kgT which gives 5 = (aj m Aj/R 4 ) with the thermal de 
Broglie wavelength given by At = (ft/ yfrnksT). 

Now for a given experimental system consisting of N - 
250000 %1 Rb atoms, we have m( 87 Rb) = 1.44 kg xlO -25 kg, 
(Oho/(2n) = 50 Hz, and an s-wave scattering length of a ^ 
5.82 x 10~ 9 m. For typical temperatures of T = 80 nK, we es- 
timate Ajja] w ~ 0.173 <K R/ai to ~ 6.37 for a condensate frac- 
tion Nq/N ~ 0.731. Thus, for typical parameter regimes, 6 2 ~ 
e 2 <s: 1 and the condensate varies on length scales much larger 
than the length scale associated with our excitations. We can, 
therefore, seek solutions in the form u\ — exp [(<p e + § e + (f>e)/d] 
and similarly for v, by using a conventional quadratic WKB 
for the excitations. Combining this with a cubic WKB ap- 
proximation for the condensate allows us to correctly resolve 
the nonlinear turning point where the thermal cloud density 
attains its maximum value. 

Motivated by the need to determine the equilibrium thermo- 
dynamic properties of a c-field simulation 11711 . we apply our 
theory to the BdG system with a Rayleigh- Jeans (RJ) distri- 
bution for nj. In the semiclassical approximation, the thermal 
population with a RJ distribution can be computed exactly as 

~ Wx) g(k,x)k 2 dk 



\nk B T 



nj 



(2n) i fia)h 
-h~(k c ,x) arctan 

+h~(k m ,x) arctan 



)8(k,x) 2 

*c00 



-/(X) 2 



h~{k c ,x) 

fcm(x) 

br (k m , x) 



- h (k c ,x) arctan 



+ h + (k m ,x) arctan 



2(k c (x)-k m (x)) 
k c (x) \ 



h + (k c ,x)l 
k m (x) ' 
h + (k m ,x) i 



,(7) 



4 



where g(k, x) = k 2 /2 + V ext (x) + 2jn (x) -fi, /(x) = yn (x), and 
/* + (A:,x) = [2g(k,x) + 2f(x)]^ 2 , h~(k,x) = [2^,x)-2/(x)] 1 / 2 . 
Therefore, we have reduced the solution of the BdG system to 
a set of algebraic equations to be solved self-consistently. We 
note that the integral in Eq. (|7]i diverges unless k c (x) is finite 
and is the well known ultraviolet catastrophe that arises in a 
classical spectrum with equipartition of energy. The choice of 
k c (x) is the single most important parameter that must be set 
in a c-fields simulation. We have set k c (x) to correspond to 
values of nj of order unity. For nj <K 1, the semi-classical ap- 
proximation breaks down and an RJ distribution significantly 
deviates from the Bose-Einstein (BE) distribution. With this 
value of «7, we have computed two cutoffs for k c . The first 
(cutoff 1) was determined using the BdG exp ression for the 
energy in the semiclassical approximation II 1611 with E c = 1 .42. 
A second cut-off using the same value of E c but with k c (x) 
now computed from the expression of the energy for a single- 
particle in a harmonic trap was also used. This choice is mo- 
tivated by the cut-offs typically used in c-field methods where 
the truncation is in terms of the single-particle basis func- 
tions. The values of k,„(x) in Eq. © are determined from 
the BdG expression for the energy by specifying a minimum 
energy E m = 0.01 as a lower bound for the integral. In Fig. QJ, 
we present numerical results of the condensate and the ther- 
mal cloud densities obtained from a self-consistent solution of 
the BdG equations with a BE and a RJ distribution using the 
method described in [18], together with the analytical results 
obtained using the first order expression for the condensate 
density. We have compared our analytical solutions obtained 
using these cut-offs against numerical results computed for the 
BdG system using generalized Laguerre basis functions. Such 
a single-particle basis provides us with more control of how to 
truncate the spectrum in our computations and has become the 
hallmark of many c-field methods |@]. In both cases, we see 
that the analytical solutions are in excellent agreement with 
the numerical results. As expected, the thermal cloud attains 
a maximum near the nonlinear turning point but since our ap- 
proximation for the condensate is smooth, no cusp arises in 
the distribution of the thermal cloud. 

Our final system is concerned with that of a nonequilibrium 
condensate. In recent years, there has been a surge of interest 
in modeling nonequilibrium systems such as exciton-polariton 
condensates flj-Sll . Recent work has shown that the steady state 
condensate density of these systems can be markedly different 
from the equilibrium Bose gases. In particular, we consider a 
2D exciton-polariton condensate under conditions studied in 
lioll . In that work, it was shown that for a given range of 
parameters, the 2D condensate profile remains radially sym- 
metric in a harmonic trap. However, the presence of pumping 
and dissipation means that a radial flow is established from 
the edges of the condensate towards the center. The flow sig- 
nificantly modifies the profile of the condensate from that ob- 
tained with a TF approximation. Here, we show how the cubic 
WKB method can be used to approximate the condensate den- 
sity in this nonequilibrium system. 
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FIG. 1: Numerical and analytical results (r measured in units of a/, ): 
(a) Numerical and analytical solutions of BE condensate at T=0; left 
and right insets show healing layer, and profiles for <p'q,<p[ respec- 
tively; (b) Condensate/thermal cloud densities at T=80nK obtained 
from numerical and analytical solutions of BdG equations with BE 
and RJ distributions; (c) Leading order approximation for density of 
2D exciton-polariton condensate; inset showing velocity within the 
condensate. 
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We begin by considering >S = atkOho and T) = crg\i/f\ 2 /fuoho, 
as in II 1 Oil - and seek a radially symmetric solution. Non- 
dimensionalising the cGPE equation in a similar way to Eq. 
([TJ, we see that a WKB ansatz applies provided a/„, <K R. 
Proceeding as above in the case of the vortex solution, we 
seek a radially symmetric solution where ifr(r) = U(x) = 
exp{[ip{x) + id(x)]/e) with r — e x to obtain two equations cor- 
responding to the real and imaginary parts 



fie 







2x = -e^Qj^-^)^ 
2 2 2 



(9) 



respectively. We now proceed as before by differentiating Eq. 
(O and using it to simplify the resulting expression. This gives 

iff [if' 2 + e<p' - ff 2 + (2^e 2v - e 4x )} + e(ffd" + e 4x - ff 2 ) 
= e 2 (<p"'/2-<p"). 

In order to approximate this equation, let us consider the re- 
gion near the origin where the exponential terms can be ne- 
glected. Moreover, the velocity at the origin must vanish since 
ff (x) = ff(r)r. We, can therefore rewrite the above equation in 
the vicinity of the origin in terms of <p'(r) as 



_ /3 eif' 2 Jiff" Iff' 



2r 2 



Now we want to enforce the boundary condition dip/dr = 
at the origin. However, from the above equation, we see that 
by neglecting all terms on the right-hand side, we can only 
enforce <p'(Q) = 0. In order to enforce dip/dr — 0, we must 
also have that <jo"(0) = 0. We can achieve this by retaining 
the last term proportional to ip'(r) appearing on the right hand 
side. Our full approximate form of <p'(x) is therefore given by 



tp ,J - + eif! - ff A + \2fLt 
+e(6'6" +e 4 - T -6>' 2 ) = 0. 



2x 



+ — 

2 



We can now solve for tp' given ff which has the solution 



(fix) = {2/ e) 



C (a-cre 2iple )t 2vle+2xf dx' 

J — CO 



-2<p/e 



(10) 



This provides an expression for the flow which is a func- 
tion of <p(x). Hence, these equations provide solutions to 
our nonequilibrium system which are expressed in terms of 
a quadrature for the density e 2 40 dx and the velocity ff (x) = 
ff(r) and must be solved self-consistently. In practice, the 
above integrals are more easily evaluated by transforming 
back to r space. 

As before, the solution we have obtained specifies the 
wavefunction upto some normalization coefficient which we 
will denote by C. However, in the case of a nonequilibrium 
condensate, we are dealing with an open system with a non- 
conserved number of particles. Therefore, in order to compu te 



the chemical potential we make use of the equation for ff (x) 
and choose C such that 



f 

xj — O 



(a- C 2 <rz 2,f ' le )e 2lfle+2x ' dx' = 0. 



(ID 



This equation defines C which once computed can be used 

to evaluate the chemical potential with the aid of Eq. (8). By 
evaluating the terms in the equation at the origin, we find ft = 

yc 2 . 

The solution of this coupled system obtained using only the 
leading order approximation for the condensate density with 
y = 0.5, a = 2.20(r pump - r), and <x = 0.15, where r pump = 4 
is the size of the pumping spot and is the step function is 
shown in Fig. lc. This corresponds to a chemical potential 
of fi = 10.97 which agrees closely with the value of 11.18 
obtained from a full numerical solution of the cGPE equa- 
tion. This numerical solution of the radial cGPE computed 
with a pseudospectral method using Laguerre-basis functions 
is shown in Fig. lc. As can be seen, the results are in excel- 
lent agreement even for this system in which the condensate 
density is significantly modified by the internal flow. 
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